
######################################################################################
#TITLE:  Introducing the Comparative Study of Electoral Systems in Tunisia: Populist Attitudes, 
#Political Preferences, and Voting Behavior
########################################################################################

#Authors:
##Ameni Mehrez , Levente Littvay , Youssef Meddeb , Bojan Todosijevic , Carsten Schneider   

##For any questions regarding this R script, please contact: Mehrez_Ameni@phd.ceu.edu



#######################################
###DATASET: CSES TUNISIA WAVE 5 #######
#######################################
##FUll codebook can be found here: https://cses.org/data-download/cses-module-5-2016-2021/

##Variable description:
https://cses.org/wp-content/uploads/2023/07/cses5_codebook_part2_variables.txt

##Parties and Leaders description:
https://cses.org/wp-content/uploads/2023/07/cses5_codebook_part3_parties_and_leaders.txt

#CSES Codebook      Election Study Code/Category
#|----------------------------------------------------------------
#  |    788001.     Abdelfattah Mourou (Ennahda Movement)
#|    788002.     Nabil Karoui (Heart of Tunisia)
#|    788003.     Abir Moussi (Free Destourian Party)
#|    788004.     Mohamed Abbou (Democratic Current)
#|    788005.     Seifeddine Makhlouf (Dignity Coalition)
#|    788007.     Youssef Chahed (Long Live Tunisia)
#|    788008.     Lotfi Mraihi (Republican People's Union)
#         |    788010.     Mehdi Jomaa (Tunisian Alternative)
#         |    788013.     Mongi Rahoui (Popular Front)
#        |    788016.     Hechmi Hamdi (Current of Love)
#       |    788019.     Moncef Marzouki (Movement Party)
#      |    788020.     Kais Saied (Independent)
#     |    788021.     Abdelkrim Zbidi (Independent)
#    |    788022.     Safi Said (Independent)
#   |    788023.     Hamma Hammami (Independent)
#  |    788024.     Elyes Fakhfakh (Democratic Forum for Labour 
# |                and Liberties)

##E2002 = Gender
##E2003 = Education
##E2010 = Income
##E2014 = Religious service
##E3004_1-7: Populism subscale
##E3005_1-5: Attitudes towards the outgroup
##E3006_1-4: Nativism subscale

library(jtools)
library(readr)
library(foreign)
library(ggplot2)
library(descr)
library(car)
library(dplyr)
library(psych)
library(effects)
library(stargazer)
library(scales)
library(sjPlot) 
library(tidyr)
library(base)
library(psych)
library(tidyr)
library(haven)
library(desc)
library(patchwork)


####DOwnload DATASET from CSES website #################
#Download link: 
https://cses.org/data-download/cses-module-5-2016-2021/
  
  
##########################################################
#############################################################

#Load data from working directory
##I am using .dta file from the CSES website

cses <- read_dta("cses5.dta")

##########################################################
#############################################################

##Get country Tunisia
TN <- subset(cses, E1006_NAM == "Tunisia")

##Get populism subscales
elite <- as.data.frame(TN[,c("E3004_1","E3004_2","E3004_3", "E3004_4","E3004_5","E3004_6","E3004_7")])
out_group <- as.data.frame(TN[,c("E3005_1","E3005_2","E3005_3","E3005_4","E3005_5")])
national <- as.data.frame(TN[,c("E3006_1", "E3006_2", "E3006_3","E3006_4")])

## sociodemographic variables
socio <- as.data.frame(TN[,c("E2002","E2003","E2010","E2014")])

##combine variables
populism <- cbind(elite, out_group, national)


##Dislike/Like parties and leaders variables
parties <- as.data.frame(TN[,c("E3017_A","E3017_B","E3017_C","E3017_D","E3017_E","E3017_F")])
leaders <- as.data.frame(TN[,c("E3018_A","E3018_B","E3018_C","E3018_D","E3018_E","E3018_F")])

##Cast a vote in the presidential elections round 1 
cast1 <- as.data.frame(TN[,c("E3012_PR_1")])

##Vote Choice: E3013_PR_1: First round presidential elections
vote1 <- as.data.frame(TN[,c("E3013_PR_1")])

###########cleaning up and recoding variables#####
##Replace DK with NAs
populism[populism[,1:16] == 7 ] <- NA
populism[populism[,1:16] == 8 ] <- NA

table(socio1$E2002)

demo <- as.data.frame(socio[,3:4])
demo[demo[,1:2] == 7 ] <- NA
demo[demo[,1:2] == 8 ] <- NA
demo[demo[,1:2] == 9 ] <- NA


socio1 <- as.data.frame(socio[,1:2])

socio1[socio1[,1:2] == 95 ] <- NA
socio1[socio1[,1:2] == 96 ] <- NA
socio1[socio1[,1:2] == 97 ] <- NA
socio1[socio1[,1:2] == 98 ] <- NA


parties[parties[,1:6] == 96 ] <- NA
parties[parties[,1:6] == 97 ] <- NA
parties[parties[,1:6] == 98 ] <- NA


leaders[leaders[,1:6] == 97 ] <- NA
leaders[leaders[,1:6] == 96 ] <- NA
leaders[leaders[,1:6] == 98 ] <- NA


##############################
#Factor Analysis ############
##############################

#This step was performed both in R and MPlus. Results reported are from the MPlus 
##The results are the same -In order to use all available observations in the dataset,
#the EFA was performed on MPlus (version 7.3) where the model is estimated using
#all datapoints using Maximum Likelihood estimation, unlike the functions ‘fa’ or 
#‘factanal’ in R that uses mean imputation or only performs EFA on complete observations. 
#The results from both softwares are very similar and can be provided upon request.

##Screeplot
#generate scree plot
populism1<- na.omit(populism[,1:16])
coro <- cor(populism1, use = "pairwise.complete.obs")
scree <- scree(coro, factors=FALSE)


##3 factor solution
fa1 <- factanal(populism1, factors = 3)
fa1


##Recoding some variables 
#(Original scale 1= strongly agree  and 5 = strongly disagree) --> I reverse that.
populism$E3004_1 <- car::recode(populism$E3004_1,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3004_2 <- car::recode(populism$E3004_2,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3004_4 <- car::recode(populism$E3004_4,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3004_5 <- car::recode(populism$E3004_5,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3004_6 <- car::recode(populism$E3004_6,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3004_7 <- car::recode(populism$E3004_7,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")


populism$E3005_1 <- car::recode(populism$E3005_1,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3005_2 <- car::recode(populism$E3005_2,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3005_4 <- car::recode(populism$E3005_4,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3005_5 <- car::recode(populism$E3005_5,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")

populism$E3006_1 <- car::recode(populism$E3006_1,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3006_2 <- car::recode(populism$E3006_2,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3006_3 <- car::recode(populism$E3006_3,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
populism$E3006_4 <- car::recode(populism$E3006_4,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")

##creating mean scores###
populism$eli <- rowMeans(populism[ , c(1,2,4,7)], na.rm=TRUE)
populism$out <- rowMeans(populism[ , c(11,12)], na.rm=TRUE)
populism$nat <- rowMeans(populism[ , c(13,14,15,16)], na.rm=TRUE)


##Changing names
names(leaders)[1]<-paste("RachedGhannouchi")
names(leaders)[2]<-paste("NabilKaroui")
names(leaders)[3]<-paste("AbirMoussi")
names(leaders)[4]<-paste("MohamedAbbou")
names(leaders)[5]<-paste("SeifeddineMakhlouf")
names(leaders)[6]<-paste("ZouheirMaghzaoui")

names(parties)[1]<-paste("Ennahda")
names(parties)[2]<-paste("QalbTounes")
names(parties)[3]<-paste("FreeDestourian")
names(parties)[4]<-paste("DemocraticCurrent")
names(parties)[5]<-paste("DignityCoalition")
names(parties)[6]<-paste("PeopleMovement")


##extracting only the 3 populist dimensions
pop_attit <- as.data.frame(populism[,17:19])


##combining data
full <- cbind(pop_attit, parties, leaders, demo, socio1,cast1, vote1)


##Changing names
names(full)[1]<-paste("Populism")
names(full)[2]<-paste("Outgroup")
names(full)[3]<-paste("Nativism")
names(full)[16]<-paste("Income")
names(full)[17]<-paste("Religiosity")

names(full)[18]<-paste("Gender")
names(full)[19]<-paste("Education")

names(full)[20]<-paste("cast1")
names(full)[21]<-paste("vote1")


###############################OLS REGRESSIONS ##########
###########Parties DV ###################################
########################################################

partya1 <- lm(Ennahda ~ Populism + Outgroup + Nativism + 
                Gender + Education + Income + 
                Religiosity, data = full)
summary(partya1)

partyb1 <- lm(QalbTounes ~ Populism + Outgroup + Nativism + 
                Gender + Education + Income + 
                Religiosity, data = full)
summary(partyb1)

partyc1 <- lm(FreeDestourian ~ Populism + Outgroup + Nativism + 
                Gender + Education + Income + 
                Religiosity, data = full)
summary(partyc1)

partyd1 <- lm(DemocraticCurrent ~ Populism + Outgroup + Nativism +
                Gender + Education + Income +
                Religiosity, data = full)
summary(partyd1)

partye1 <- lm(DignityCoalition ~ Populism + Outgroup + Nativism +
                Gender + Education + Income + 
                Religiosity, data = full)
summary(partye1)

partyf1 <- lm(PeopleMovement ~ Populism + Outgroup + Nativism +
                Gender + Education + Income + 
                Religiosity, data = full)
summary(partyf1)


##Plot results (without democracy variable)
aa1 <- plot_summs(partya1, partyb1, partyc1, partyd1, partye1, partyf1, 
                  scale = TRUE, robust = TRUE,
                  point.size = 3, #inner_ci_level = .3,
                  point.shape = T,
                  model.names = c("Ennahda", "Qalb Tounes", 
                                  "Free Destourian", "Democratic Current",
                                  "Dignity Coalition", "People Movement"))
aa1

ff <- aa1 + xlab("Coefficient estimate") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.line.x = element_line(colour = "black", linetype = "solid"),
        axis.ticks.x = element_line(colour = "black", linetype = "solid"),
        axis.text=element_text(size=8),
        axis.title=element_text(size=12,face="bold"))

ff

ggsave(plot = ff ,"fig1.tiff", dpi=600, units = "px", 
       height=3862, width=3700)



ggsave(plot = ff ,"fig1.pdf", dpi=600, units = "px", 
       height=3862, width=3700)

###create regression table in latex
stargazer(partya1, partyb1, partyc1, 
          partyd1, partye1, partyf1, type = "latex",
          omit.stat=c("LL","ser","f"))


###################################################################
#####Leaders as DV#################################################

leadera <- lm(RachedGhannouchi ~ Populism + Outgroup + Nativism + 
                Gender + Education + Income + Religiosity, data = full)
summary(leadera)

leaderb <- lm(NabilKaroui ~ Populism + Outgroup + Nativism + 
                Gender + Education + Income + Religiosity, data = full)
summary(leaderb)

leaderc <- lm(AbirMoussi ~ Populism + Outgroup + Nativism + 
                Gender + Education + Income + Religiosity, data = full)
summary(leaderc)

leaderd <- lm(MohamedAbbou ~ Populism + Outgroup + Nativism +
                Gender + Education + Income + Religiosity, data = full)
summary(leaderd)

leadere <- lm(SeifeddineMakhlouf ~ Populism + Outgroup + Nativism +
                Gender + Education + Income + Religiosity, data = full)
summary(leadere)

leaderf <- lm(ZouheirMaghzaoui ~ Populism + Outgroup + Nativism + 
                Gender + Education + Income + Religiosity, data = full)
summary(leaderf)


#Plot models
bb <- plot_summs(leadera, leaderb, leaderc, leaderd, 
                 leadere, leaderf, point.size = 3,
                 point.shape = T,
                 scale = TRUE, robust = TRUE,
                 model.names = c("Rached Ghannouchi", 
                                 "Nabil Karoui", 
                                 "Abir Moussi", 
                                 "Mohamed Abbou",
                                 "Seifeddine Makhlouf",
                                 "Zouheir Maghzaoui"))

bb


hh <- bb + xlab("Coefficient estimate") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.line.x = element_line(colour = "black", linetype = "solid"),
        axis.ticks.x = element_line(colour = "black", linetype = "solid"),
        axis.text=element_text(size=8),
        axis.title=element_text(size=12,face="bold"))

hh

##Save plot for leaders ###
ggsave(plot = hh ,"fig2.tiff", dpi=600,units = "px", 
       height=3362, width=3700)

ggsave(plot = hh ,"fig2.pdf", dpi=600,units = "px", 
       height=3362, width=3700)

stargazer(leadera, leaderb, leaderc, 
          leaderd, leadere, leaderf, type = "latex",
          omit.stat=c("LL","ser","f"))



##############################
##Recode  presidential votes
############################

table(full$cast1)
full$cast1 <- car::recode(full$cast1,"93=NA; 98=NA")
table(full$vote1)

#Vote choice
full$vote1 <- car::recode(full$vote1,"788001=0; 
                               788002=0; 788003= 0; 788004 = 0;
                               788005 = 0; 788007 =0; 788008 =0;
                               788010 =0; 788013 =0; 788016=0;
                               788019 = 0; 788021=0; 788022=0;
                               788023 =0; 788024 =0; 999992=0;
                               788020=1;
                               999993=NA;999997=NA;
                               999998=NA;999999=NA")


##Remove double labels
fulls <- sapply(full, haven::zap_labels)
fulls <- as.data.frame(fulls)

#####################################################
##Logit models########################################
#######################################################

cast_vote <- glm(cast1 ~ Populism + Outgroup + Nativism +
                   Gender + Education + Income + Religiosity, 
                 data = fulls, family = "binomial")
summary(cast_vote)
plot_coefs(cast_vote)
exp(coef(cast_vote))


#################Logit model ###############
########Voting for Kais vs. others
presi_vote1 <- glm(vote1 ~ Populism + Outgroup + Nativism +
                     Gender + Education + Income + Religiosity, 
                   data = fulls, family = "binomial")
summary(presi_vote1)

exp(coef(presi_vote1))

1.36-1 #36%

##Save into stargazer table
stargazer(cast_vote, presi_vote1, type = "latex")

##Plot exponents
cc <- plot_summs(cast_vote, presi_vote1, exp = T, point.size = 3,
                 point.shape = T,
                 robust = TRUE,
                 model.names = c("Voters", "Kais Saied"))
cc

cc1 <- cc + xlab("Odds ratios") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.line.x = element_line(colour = "black", linetype = "solid"),
        axis.ticks.x = element_line(colour = "black", linetype = "solid"),
        axis.text=element_text(size=8),
        axis.title=element_text(size=12,face="bold"))
cc1



ggsave(plot = cc1 ,"fig3.tiff", dpi=600,units = "px", 
       height=3362, width=3700)

ggsave(plot = cc1 ,"fig3.pdf", dpi=600,units = "px", 
       height=3362, width=3700)

###Plot descriptive stats - TUNISIA
##Politicians do not care about the people

pop2 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(210,88,77,184,888)
)


pop2 <- pop2 %>% mutate(freq = scales::label_percent()(count/sum(count)))

po2 <- ggplot(pop2, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())
#axis.text.y = element_text("count"))
#axis.text.y = element_text(size = 6)) 

po2

po22 <- po2 + ggtitle("Most politicians do not care about the people")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po22
##what people call compromise is really just selling out on one's principles
table(populism$E3004_1)

pop1 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(390,175,251,261,269)
)

pop1 <- pop1 %>% mutate(freq = scales::label_percent()(count/sum(count)))

#pop1$data$name <- factor(pop2$data$name, levels = c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po1 <- ggplot(pop1, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())
#axis.text.y = element_text("count"))
#axis.text.y = element_text(size = 6)) 

po1
po11 <- po1 + ggtitle("What people call compromise is really just selling out on one's \n principles")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po11
##pop3 most politicians are trustworthy


pop3 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(897,198,86,181,72)
)
pop3

pop3 <- pop3 %>% mutate(freq = scales::label_percent()(count/sum(count)))

po3 <- ggplot(pop3, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

po3
po33 <- po3 + ggtitle("Most politicians are trustworthy")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po33



###Politicians are the main problem in Tunisia

pop4 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(183,130,108,188,813)
)
pop4

pop4 <- pop4 %>% mutate(freq = scales::label_percent()(count/sum(count)))

po4 <- ggplot(pop4, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

po4
po44 <- po4 + ggtitle("Politicians are the main problem in Tunisia")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po44


###Having a strong leader in government is good for Tunisia even if the leader bends the rules to get things done


pop5 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(308,102,106,263,662)
)


pop5 <- pop5 %>% mutate(freq = scales::label_percent()(count/sum(count)))
pop5

po5 <- ggplot(pop5, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

po5
po55 <- po5 + ggtitle("Having a strong leader in government is good for Tunisia even if the leader \n bends the rules to get things done")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po55

###po6 The people and not the politicians should make our most important policy decisions

pop6 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(204,154,131,278,677)
)


pop6 <- pop6 %>% mutate(freq = scales::label_percent()(count/sum(count)))
pop6

po6 <- ggplot(pop6, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

po6
po66 <- po6 + ggtitle("The people and not the politicians should make our most \n important policy decisions")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po66


### Most politicians care only about the interests of the rich and powerful
pop7 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(199,91,78,189,894)
)


pop7 <- pop7 %>% mutate(freq = scales::label_percent()(count/sum(count)))
pop7

po7 <- ggplot(pop7, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

po7
po77 <- po7 + ggtitle("Most politicians care only about the interests of the rich \n and powerful")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

po77

##################################################

po11 + po22 + po33 + po44 + po55 + po66 + po77 + plot_layout(ncol=2, nrow = 4)

ggsave("pag1.png", dpi=500,units = "px", 
       height=8862, width=6200)



######Move to attitudes towards the outgroup ####
### Minorities should adapt to the customs and traditions of Tunisia
table(populism$E3005_1)
atti1 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(193,125,179,350,580))

atti1 <- atti1 %>% mutate(freq = scales::label_percent()(count/sum(count)))
atti1

att1 <- ggplot(atti1, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

att1
att11 <- att1 + ggtitle("Minorities should adapt to the customs and traditions of Tunisia")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

att11



#The will of the majority should always prevail, even over the rights of minorities
table(populism$E3005_2)
atti2 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(484,182,159,254,342))

atti2 <- atti2 %>% mutate(freq = scales::label_percent()(count/sum(count)))
atti2

att2 <- ggplot(atti2, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

att2
att22 <- att2 + ggtitle("The will of the majority should always prevail, even over the rights \n of minorities")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

att22



#Immigrants are generally good for Tunisia's economy
table(populism$E3005_3)

atti3 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(246,116,94,349,635))

atti3 <- atti3 %>% mutate(freq = scales::label_percent()(count/sum(count)))
atti3

att3 <- ggplot(atti3, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

att3
att33 <- att3 + ggtitle("Immigrants are generally good for Tunisia's economy")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

att33


## Tunisia's culture is generally harmed by immigrants
table(populism$E3005_4)

atti4 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(517,294,161,210,245))

atti4 <- atti4 %>% mutate(freq = scales::label_percent()(count/sum(count)))
atti4

att4 <- ggplot(atti4, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

att4
att44 <- att4 + ggtitle("Tunisia's culture is generally harmed by immigrants")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

att44

#Immigrants increase crime rates in Tunisia
table(populism$E3005_5)

atti5 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(616,267,171,155,203))

atti5 <- atti5 %>% mutate(freq = scales::label_percent()(count/sum(count)))
atti5

att5 <- ggplot(atti5, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

att5
att55 <- att5 + ggtitle("Immigrants increase crime rates in Tunisia")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

att55

att11 + att22 + att33 + att44 + att55 + plot_layout(ncol=2, nrow = 3)

ggsave("pag2.png", dpi=500,units = "px", 
       height=8862, width=6200)


#####################################################################################
####nativism scale: Important for being truly Tunisian: To have been born in Tunisia
table(cses$E3006_1)

table(populism$E3006_1)

nati1 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(0,145,185,162,973))

nati1 <- nati1 %>% mutate(freq = scales::label_percent()(count/sum(count)))
nati1

nat1 <- ggplot(nati1, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

nat1
nat11 <- nat1 + ggtitle("Important for being truly Tunisian: To have been born in Tunisia")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

nat11


#### To have Tunisian ancestry
table(cses$E3006_2)

table(populism$E3006_2)

nati2 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(0,75,122,171,1102))

nati2 <- nati2 %>% mutate(freq = scales::label_percent()(count/sum(count)))
nati2

nat2 <- ggplot(nati2, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

nat2
nat22 <- nat2 + ggtitle("Important for being truly Tunisian: To have Tunisian ancestry")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

nat22

#### To be able to speak Tunisian
table(cses$E3006_3)

table(populism$E3006_3)

nati3 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(0,98,203,217,947))

nati3 <- nati3 %>% mutate(freq = scales::label_percent()(count/sum(count)))
nati3

nat3 <- ggplot(nati3, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

nat3
nat33 <- nat3 + ggtitle("Important for being truly Tunisian: To be able to speak Tunisian")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

nat33


#### To follow Tunisia's customs and traditions
table(cses$E3006_4)

nati4 <- data.frame(
  name=c("Strongly disagree","Disagree","Neither A nor D","Agree","Strongly agree") ,  
  count=c(0,77,126,242,1024))

nati4 <- nati4 %>% mutate(freq = scales::label_percent()(count/sum(count)))
nati4

nat4 <- ggplot(nati4, aes(x= name, y=count)) + 
  geom_bar(stat="identity") +
  geom_text(aes(label=freq), hjust=0.5, vjust= -0.4 , size=3)+
  #ylab("count")+
  #coord_flip() +
  theme_bw()+ theme(axis.text = element_text(size = 8),
                    axis.title.x = element_blank(),
                    axis.text.y =  element_blank(),
                    axis.ticks.y =  element_blank(),
                    axis.title.y = element_blank())

nat4
nat44 <- nat4 + ggtitle("Important for being truly Tunisian: To follow Tunisia's customs and traditions")+ 
  scale_x_discrete(limits=c("Strongly agree","Agree", "Neither A nor D","Disagree","Strongly disagree"))

nat44

nat11 + nat22 + nat33 + nat44 + plot_layout(ncol=2, nrow = 2)

ggsave("pag3.png", dpi=500,units = "px", 
       height=8362, width=6300)





###################################################
###comparison across countries###################
################################################

##Subset some countries only from the whole dataset 
cses2 <- subset(cses, E1006_NAM == "Tunisia" | E1006_NAM == "France" |
                  E1006_NAM == "United States of America"| E1006_NAM == "Brazil" |
                  E1006_NAM == "Great Britain"| E1006_NAM == "Hungary"| E1006_NAM == "Turkey"|
                  E1006_NAM == "Israel")


cses3 <- as.data.frame(cses2[,92:118])
coun <- as.data.frame(cses2[,c("E1006_NAM")])


cses3[cses3[,1:27] == 7] <- NA
cses3[cses3[,1:27] == 8] <- NA
cses3[cses3[,1:27] == 9] <- NA

cses3 <- cbind(coun, cses3)


##Recode key populism items
cses3$E3004_1 <- car::recode(cses3$E3004_1,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3004_2 <- car::recode(cses3$E3004_2,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3004_4 <- car::recode(cses3$E3004_4,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3004_5 <- car::recode(cses3$E3004_5,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3004_6 <- car::recode(cses3$E3004_6,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3004_7 <- car::recode(cses3$E3004_7,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")


cses3$E3005_1 <- car::recode(cses3$E3005_1,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3005_2 <- car::recode(cses3$E3005_2,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3005_4 <- car::recode(cses3$E3005_4,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3005_5 <- car::recode(cses3$E3005_5,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")

cses3$E3006_1 <- car::recode(cses3$E3006_1,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3006_2 <- car::recode(cses3$E3006_2,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3006_3 <- car::recode(cses3$E3006_3,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")
cses3$E3006_4 <- car::recode(cses3$E3006_4,"1=5;2=4;3=3;4=2;5=1;7=NA;8=NA;9=NA")


means1 <- cses3 %>% 
  group_by(E1006_NAM) %>% 
  summarise_at(vars("E3004_1", "E3004_2", "E3004_3", "E3004_4", "E3004_5", "E3004_6", "E3004_7"), mean,na.rm = T)

#summarise(E3004_1 =mean(E3004_1, na.rm=T)) 
means1

means1$elite <- rowMeans(means1[ , c("E3004_1", "E3004_2", "E3004_3", "E3004_4", "E3004_5", "E3004_6", "E3004_7")], na.rm=TRUE)


# View the model summary
summary(model)

means1$E1006_NAM <- as.factor(means1$E1006_NAM)

means1$E1006_NAM <- factor(means1$E1006_NAM, levels = c("Brazil", "Tunisia", "Hungary", "Turkey", "France", "United States of America","Israel", "Great Britain" ))


table(means1$E1006_NAM)
means <- as.data.frame(cses3[,1:9])
means <- as.data.frame(means[,-3])


# RUN THIS CODE if the next script does not work: detach(package:plyr)
#library(dplyr)
hh1 <- means %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean1 = mean(E3004_1, na.rm = T),
            lci = t.test(E3004_1, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3004_1, conf.level = 0.95)$conf.int[2])

hh2 <- means %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean2 = mean(E3004_2, na.rm = T),
            lci = t.test(E3004_2, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3004_2, conf.level = 0.95)$conf.int[2])

hh3 <- means %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean3 = mean(E3004_3, na.rm = T),
            lci = t.test(E3004_3, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3004_3, conf.level = 0.95)$conf.int[2])



hh4 <- means %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean4 = mean(E3004_4, na.rm = T),
            lci = t.test(E3004_4, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3004_4, conf.level = 0.95)$conf.int[2])


hh5 <- means %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean5 = mean(E3004_5, na.rm = T),
            lci = t.test(E3004_5, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3004_5, conf.level = 0.95)$conf.int[2])


hh6 <- means %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean6 = mean(E3004_6, na.rm = T),
            lci = t.test(E3004_6, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3004_6, conf.level = 0.95)$conf.int[2])


hh7 <- means %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean7 = mean(E3004_7, na.rm = T),
            lci = t.test(E3004_7, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3004_7, conf.level = 0.95)$conf.int[2])


hh <- cbind(hh1,hh2,hh3,hh4,hh5,hh6,hh7)

hh$elite <- rowMeans(hh[ , c("mean1", "mean2", "mean3", "mean4", "mean5", "mean6", "mean7")], na.rm=TRUE)

hh$lower <- rowMeans(hh[ , c(3,7,11,15,19,23,27)])

hh$upper <- rowMeans(hh[ , c(4,8,12,16,20,24,28)])


final1 <- as.data.frame(hh[,c("E1006_NAM","elite","lower", "upper")])

final1$E1006_NAM <- as.factor(final1$E1006_NAM)

final1$E1006_NAM <- factor(final1$E1006_NAM, levels = c("Brazil", "Tunisia", "Hungary", "Turkey", "France", "United States of America","Israel", "Great Britain" ))

p <- ggplot(final1, aes(x = E1006_NAM, y = elite))+
  geom_point(size = 2, color = "red") +
  geom_errorbar(aes(ymin=lower, ymax=upper),width = 0.1, size = 0.8)+ 
  labs(x = "Countries", y = "Populist Attitudes") +
  theme_bw()+ theme(axis.text = element_text(size = 6),
                    axis.title.y = element_blank(),
                    axis.text.y = element_text(size = 8))

pp<- p + coord_flip()
pp

ggsave(plot = pp ,"countries1.png", dpi=500,units = "px", 
       height=2362, width=3500)


###Attitudes towards the out-group

means2 <- as.data.frame(cses3[,c(1,11,12,13,14,15)])


kk1 <- means2 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean1 = mean(E3005_1, na.rm = T),
            lci = t.test(E3005_1, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3005_1, conf.level = 0.95)$conf.int[2])

kk2 <- means2 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean2 = mean(E3005_2, na.rm = T),
            lci = t.test(E3005_2, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3005_2, conf.level = 0.95)$conf.int[2])

kk3 <- means2 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean3 = mean(E3005_3, na.rm = T),
            lci = t.test(E3005_3, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3005_3, conf.level = 0.95)$conf.int[2])



kk4 <- means2 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean4 = mean(E3005_4, na.rm = T),
            lci = t.test(E3005_4, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3005_4, conf.level = 0.95)$conf.int[2])


kk5 <- means2 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean5 = mean(E3005_5, na.rm = T),
            lci = t.test(E3005_5, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3005_5, conf.level = 0.95)$conf.int[2])




kk <- cbind(kk1,kk2,kk3,kk4,kk5)

kk$outgroup <- rowMeans(kk[ , c("mean1", "mean2", "mean3", "mean4", "mean5")], na.rm=TRUE)

kk$lower <- rowMeans(kk[ , c(3,7,11,15,19)])

kk$upper <- rowMeans(kk[ , c(4,8,12,16,20)])

final2 <- as.data.frame(kk[,c("E1006_NAM","outgroup","lower", "upper")])

final2$E1006_NAM <- as.factor(final2$E1006_NAM)
final2$E1006_NAM <- factor(final2$E1006_NAM, levels = c("Hungary","Turkey", "France","Great Britain","Israel", "Brazil", "Tunisia","United States of America"))


q <- ggplot(final2, aes(x = E1006_NAM, y = outgroup))+
  geom_point(size = 2, color = "red") +
  geom_errorbar(aes(ymin=lower, ymax=upper),width = 0.1, size = 0.8)+ 
  labs(x = "Countries", y = "Attitudes towards the outgroup") +
  theme_bw()+ theme(axis.text = element_text(size = 6),
                    axis.title.y = element_blank(),
                    axis.text.y = element_text(size = 8))

qq<- q + coord_flip()
qq

ggsave(plot = qq ,"countries2.png", dpi=500,units = "px", 
       height=2362, width=3500)


####Nativism scale ###########

###Attitudes towards the out-group

means3 <- as.data.frame(cses3[,c(1,16,17,18,19)])


jj1 <- means3 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean1 = mean(E3006_1, na.rm = T),
            lci = t.test(E3006_1, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3006_1, conf.level = 0.95)$conf.int[2])

jj2 <- means3 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean2 = mean(E3006_2, na.rm = T),
            lci = t.test(E3006_2, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3006_2, conf.level = 0.95)$conf.int[2])

jj3 <- means3 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean3 = mean(E3006_3, na.rm = T),
            lci = t.test(E3006_3, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3006_3, conf.level = 0.95)$conf.int[2])


jj4 <- means3 %>% 
  group_by(E1006_NAM) %>% 
  summarise(mean4 = mean(E3006_4, na.rm = T),
            lci = t.test(E3006_4, conf.level = 0.95)$conf.int[1],
            uci = t.test(E3006_4, conf.level = 0.95)$conf.int[2])


jj <- cbind(jj1,jj2,jj3,jj4)

jj$nativism <- rowMeans(jj[ , c("mean1", "mean2", "mean3", "mean4")], na.rm=TRUE)

jj$lower <- rowMeans(jj[ , c(3,7,11,15)])

jj$upper <- rowMeans(jj[ , c(4,8,12,16)])

final3 <- as.data.frame(jj[,c("E1006_NAM","nativism","lower", "upper")])

final3$E1006_NAM <- as.factor(final3$E1006_NAM)
final3$E1006_NAM <- factor(final3$E1006_NAM, levels = c("Tunisia","Hungary", "Brazil","Turkey", "France","Great Britain", "Israel","United States of America"))


l <- ggplot(final3, aes(x = E1006_NAM, y = nativism))+
  geom_point(size = 2, color = "red") +
  geom_errorbar(aes(ymin=lower, ymax=upper),width = 0.1, size = 0.8)+ 
  labs(x = "Countries", y = "Nativism") +
  theme_bw()+ theme(axis.text = element_text(size = 6),
                    axis.title.y = element_blank(),
                    axis.text.y = element_text(size = 8))

ll <- l + coord_flip()
ll

ggsave(plot = ll ,"countries3.png", dpi=500,units = "px", 
       height=2362, width=3500)


o <- pp + qq + ll + plot_layout(ncol=1, nrow = 3)+
  plot_annotation(title ="" ,
                  theme = theme(plot.title = element_text(size = 13)))

o

#######################################
###Plotting means across countries####
#####################################

library(forcats)
plot1 <-  means1 %>% 
  arrange(elite) %>%
  ungroup() %>%
  mutate(names = fct_reorder(E1006_NAM, elite), .desc = F) %>%
  ggplot(aes(x=E1006_NAM,y=elite)) + 
  geom_bar(position=position_dodge(),stat="identity", fill = "cornflowerblue") +
  coord_flip() +
  ylab("Populist Attitudes")+
  theme_bw()+ theme(axis.text = element_text(size = 6),
                    axis.title.y = element_blank(),
                    axis.text.y = element_text(size = 6)) 
plot1


##outgroup
means2 <- cses3 %>% 
  group_by(E1006_NAM) %>% 
  summarise_at(vars("E3005_1", "E3005_2", "E3005_3", "E3005_4", "E3005_5"), mean, na.rm = T)
#summarise(E3004_1 =mean(E3004_1, na.rm=T)) 
means2

means2$outgroup <- rowMeans(means2[ , c("E3005_1", "E3005_2", "E3005_3", "E3005_4", "E3005_5")], na.rm=TRUE)


means2$E1006_NAM <- as.factor(means2$E1006_NAM)

means2$E1006_NAM <- factor(means2$E1006_NAM, levels = c("Hungary","Turkey", "France","Great Britain","Israel", "Brazil", "Tunisia","United States of America"))


plot2 <-  means2 %>% 
  arrange(outgroup) %>%
  ungroup() %>%
  mutate(names = fct_reorder(E1006_NAM, outgroup), .desc = F) %>%
  ggplot(aes(x=E1006_NAM,y=outgroup)) + 
  geom_bar(position=position_dodge(),stat="identity", fill = "cornflowerblue") +
  coord_flip() +
  ylab("Attitudes towards the outgroup")+
  theme_bw()+ theme(axis.text = element_text(size = 6),
                    axis.title.y = element_blank(),
                    axis.text.y = element_text(size = 6)) 
plot2


##Nativism
means3 <- cses3 %>% 
  group_by(E1006_NAM) %>% 
  summarise_at(vars("E3006_1", "E3006_2", "E3006_3", "E3006_4"), mean, na.rm = T)
#summarise(E3004_1 =mean(E3004_1, na.rm=T)) 
means3

means3$nativism <- rowMeans(means3[ , c("E3006_1", "E3006_2", "E3006_3", "E3006_4")], na.rm=TRUE)


means3$E1006_NAM <- as.factor(means3$E1006_NAM)

means3$E1006_NAM <- factor(means3$E1006_NAM, levels = c("Tunisia","Hungary", "Brazil","Turkey", "France","Great Britain", "Israel","United States of America"))


plot3 <-  means3 %>% 
  arrange(nativism) %>%
  ungroup() %>%
  mutate(names = fct_reorder(E1006_NAM, nativism), .desc = F) %>%
  ggplot(aes(x=E1006_NAM,y=nativism)) + 
  geom_bar(position=position_dodge(),stat="identity", fill = "cornflowerblue") +
  coord_flip() +
  ylab("Nativism")+
  theme_bw()+ theme(axis.text = element_text(size = 6),
                    axis.title.y = element_blank(),
                    axis.text.y = element_text(size = 6)) 
plot3




ggsave(plot = plot1 ,"comp1.png", dpi=500,units = "px", 
       height=2362, width=3700)



ggsave(plot = plot2 ,"comp2.png", dpi=500,units = "px", 
       height=2362, width=3700)


ggsave(plot = plot3 ,"comp3.png", dpi=500,units = "px", 
       height=2362, width=3700)


################################################################
###################### Making maps Tunisian data ###############

# Load required libraries
library(maps)
library(ggplot2)
library(readr)
library(raster)
library(scales)
library(scales)


##choloplepth

dt<- read_csv("tn_gouv.csv") 

m_gouv2<- raster::getData(name="GADM",  country="TUN", level=1)
m_gouv2@data$HASC_1[1:4]
tn_gouv_fr <- fortify(m_gouv2, region = "HASC_1")

i=match(tn_gouv_fr$id,dt$HASC_1)
tn_gouv_fr$sample=dt$sample[i]
colnames(tn_gouv_fr)


p2 <-ggplot(tn_gouv_fr, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(fill=sample),color = "black")+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()


p2<-p2+theme(axis.line=element_blank(),
             axis.text.x=element_blank(),
             axis.text.y=element_blank(),
             axis.ticks=element_blank(),
             axis.title.x=element_blank(),
             axis.title.y=element_blank(),
             panel.grid.major = element_blank(),
             panel.grid.minor = element_blank(),
             panel.border = element_blank(),
             panel.background = element_blank())
p2<-p2+ theme(legend.position="right")
p3 <- p2+ scale_fill_viridis_c(labels= scales::comma)
p3

ggsave(plot = p3 ,"map1.png", dpi=500,units = "px", 
       height=2362, width=1800)



##do the same for population
m_gouv3 <- raster::getData(name="GADM",  country="TUN", level=1)
m_gouv3@data$HASC_1[1:4]
tn_gouv_fr <- fortify(m_gouv3, region = "HASC_1")

i=match(tn_gouv_fr$id,dt$HASC_1)
tn_gouv_fr$population=dt$population[i]
colnames(tn_gouv_fr)


p3 <-ggplot(tn_gouv_fr, aes(x=long, y=lat, group=group)) +
  geom_polygon(aes(fill=population),color = "black")+
  labs(x="",y="")+ theme_bw()+
  coord_fixed()

p3<-p3+theme(axis.line=element_blank(),
             axis.text.x=element_blank(),
             axis.text.y=element_blank(),
             axis.ticks=element_blank(),
             axis.title.x=element_blank(),
             axis.title.y=element_blank(),
             panel.grid.major = element_blank(),
             panel.grid.minor = element_blank(),
             panel.border = element_blank(),
             panel.background = element_blank())
p3<-p3+ theme(legend.position="right")
p4<- p3 + scale_fill_viridis_c(labels= scales::comma)
p4

ggsave(plot = p4 ,"map2.png", dpi=500,units = "px", 
       height=2362, width=1800)



p5 <- p3 + p4 + plot_layout(ncol=2, nrow = 1)+
  plot_annotation(title ="Choropleth Maps of Sample vs. Population Distributions" ,
                  theme = theme(plot.title = element_text(size = 13)))

p5

ggsave(plot = p5 ,"maptogether.png", dpi=500,units = "px", 
       height=2362, width=3600)


#################################################################################
###Do the same analysis by including only people who voted in parliament election
##subset voters and exclude non-voters among party supporters#############
parl <- as.data.frame(TN[,"E3012_LH"])
fulla <- cbind(full, parl)

fulla <- sapply(fulla, haven::zap_labels)
fulla <- as.data.frame(fulla)


fulla$E3012_LH <- car::recode(fulla$E3012_LH,"93=NA;97=NA ;98=NA; 99=NA")

table(fulla$E3012_LH)

names(fulla)[22]<-paste("parl")


fulla <- subset(fulla, fulla$parl == 1)

partya01 <- lm(Ennahda ~ Populism + Outgroup + Nativism + 
                 Gender + Education + Income + 
                 Religiosity, data = fulla)
summary(partya01)

partyb01 <- lm(QalbTounes ~ Populism + Outgroup + Nativism + 
                 Gender + Education + Income + 
                 Religiosity, data = fulla)
summary(partyb01)

partyc01 <- lm(FreeDestourian ~ Populism + Outgroup + Nativism + 
                 Gender + Education + Income + 
                 Religiosity, data = fulla)
summary(partyc01)

partyd01 <- lm(DemocraticCurrent ~ Populism + Outgroup + Nativism +
                 Gender + Education + Income +
                 Religiosity, data = fulla)
summary(partyd01)

partye01 <- lm(DignityCoalition ~ Populism + Outgroup + Nativism +
                 Gender + Education + Income + 
                 Religiosity, data = fulla)
summary(partye01)

partyf01 <- lm(PeopleMovement ~ Populism + Outgroup + Nativism +
                 Gender + Education + Income + 
                 Religiosity, data = fulla)
summary(partyf01)


##Plot results (without democracy variable)
aa2 <- plot_summs(partya01, partyb01, partyc01, partyd01, partye01, partyf01, 
                  scale = TRUE, robust = TRUE,
                  point.size = 3, #inner_ci_level = .3,
                  point.shape = T,
                  model.names = c("Ennahda", "Qalb Tounes", 
                                  "Free Destourian", "Democratic Current",
                                  "Dignity Coalition", "People Movement"))

aa2


ff2 <- aa2 + xlab("Coefficient estimate") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.line.x = element_line(colour = "black", linetype = "solid"),
        axis.ticks.x = element_line(colour = "black", linetype = "solid"),
        axis.text=element_text(size=8),
        axis.title=element_text(size=12,face="bold"))

ff2
#pp <- aa1 +            # Change margins of ggplot2 plot
# theme(plot.margin = unit(c(0.5, 1.5,0.5,1), "cm"))

ggsave(plot = ff ,"fig1.tiff", dpi=600, units = "px", 
       height=3862, width=3700)


stargazer(partya01, partyb01, partyc01, 
          partyd01, partye01, partyf01, type = "latex",
          omit.stat=c("LL","ser","f"))


##Re-run regression analysis with voters among leader supporters
table(full$presi)
fullb <- subset(fulla, full$cast1 == 1)

leadera1 <- lm(RachedGhannouchi ~ Populism + Outgroup + Nativism + 
                 Gender + Education + Income + Religiosity, data = fullb)
summary(leadera1)

leaderb1 <- lm(NabilKaroui ~ Populism + Outgroup + Nativism + 
                 Gender + Education + Income + Religiosity, data = fullb)
summary(leaderb1)

leaderc1 <- lm(AbirMoussi ~ Populism + Outgroup + Nativism + 
                 Gender + Education + Income + Religiosity, data = fullb)
summary(leaderc1)

leaderd1 <- lm(MohamedAbbou ~ Populism + Outgroup + Nativism +
                 Gender + Education + Income + Religiosity, data = fullb)
summary(leaderd1)

leadere1 <- lm(SeifeddineMakhlouf ~ Populism + Outgroup + Nativism +
                 Gender + Education + Income + Religiosity, data = fullb)
summary(leadere1)

leaderf1 <- lm(ZouheirMaghzaoui ~ Populism + Outgroup + Nativism + 
                 Gender + Education + Income + Religiosity, data = fullb)
summary(leaderf1)


bb1 <- plot_summs(leadera1, leaderb1, leaderc1, leaderd1, 
                  leadere1, leaderf1, point.size = 3,
                  point.shape = T,
                  scale = TRUE, robust = TRUE,
                  model.names = c("Rached Ghannouchi", 
                                  "Nabil Karoui", 
                                  "Abir Moussi", 
                                  "Mohamed Abbou",
                                  "Seifeddine Makhlouf",
                                  "Zouheir Maghzaoui"))

bb1


hh1 <- bb1 + xlab("Coefficient estimate") + ylab("") +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
  theme(axis.line.x = element_line(colour = "black", linetype = "solid"),
        axis.ticks.x = element_line(colour = "black", linetype = "solid"),
        axis.text=element_text(size=8),
        axis.title=element_text(size=12,face="bold"))

hh1

##Save plot for leaders ###
ggsave(plot = hh1 ,"fig2.tiff", dpi=600,units = "px", 
       height=3362, width=3700)



stargazer(leadera1, leaderb1, leaderc1, 
          leaderd1, leadere1, leaderf1, type = "latex",
          omit.stat=c("LL","ser","f"))

#####################################
###Descriptive stats
####################################

cses5 <- TN %>%
  dplyr::select("E2002","E2003","E2010","E2014", "E3001","E3002","E3003",
                "E3017_A","E3017_B","E3017_C","E3017_D","E3017_E","E3017_F",
                "E3018_A","E3018_B","E3018_C","E3018_D","E3018_E","E3018_F")

cses5[cses5[,1:19] == 95 ] <- NA
cses5[cses5[,1:19] == 96 ] <- NA
cses5[cses5[,1:19] == 97 ] <- NA
cses5[cses5[,1:19] == 98 ] <- NA
cses5[cses5[,1:19] == 99 ] <- NA


cses5_1 <- as.data.frame(cses5[,1:7])

cses5_1[cses5_1[,1:7] == 7 ] <- NA
cses5_1[cses5_1[,1:7] == 8 ] <- NA
cses5_1[cses5_1[,1:7] == 9 ] <- NA

cses5_2 <- as.data.frame(cses5[,8:19])

combined <- cbind(cses5_1, cses5_2)

##recode variables
table(combined$E3017_B)

combined$E3001 <- car::recode(combined$E3001,"1=4;2=3;3=2;4=1")
combined$E3002 <- car::recode(combined$E3002,"1=4;2=3;3=2;4=1")
combined$E3003 <- car::recode(combined$E3003,"1=5;2=4;3=3;4=2;5=1")


dat <- as.data.frame(describe(combined[,1:19]))

write.csv(dat, "descriptive.csv")
